Adaptive weighted stochastic gradient descent

ABSTRACT

A system and method for solving an optimization problem are disclosed which are based on adaptive weighted stochastic gradient descent. The method includes receiving data to be sampled in optimizing an integral of a second function and, for a number of iterations, sampling the data according to a modified distribution function which depends on at least one parameter, updating at least one parameter of the second function, and updating at least one parameter of the modified distribution function. In a subsequent one of the iterations, the modified distribution function used in sampling the data is based on the updated at least one parameter of the distribution function. Information is output, which can be based on the integral function with the at least one updated parameter.

BACKGROUND

The exemplary embodiment relates to optimization of integral functionsand finds particular application in connection with a system and methodwhich used an adapted weighted stochastic gradient descent approach toupdate a probability distribution function used for sampling data duringthe optimization.

Many problems involve optimizing a function (e.g., minimizing a costfunction or maximizing a future reward) with respect to some parameters(such as the investment level on each different asset). The function tooptimize is often an integral (for example, an expectation, i.e., theaverage over multiple possible outcomes). When this integral is nottractable, i.e., when it is hard to compute, the optimization problem isitself difficult. In many real-world problems, such as those involvingintractable integrals, such as the averaging on combinatorial spaces,non-Gaussian integrals, etc., the solution can be expressed as anintegral of a function of a sampled value and a gradient of thefunction. Stochastic Gradient Descent (SGD) is an optimization techniquethat can optimize an integral without having to compute it. However, SGDis known to be slow to converge.

Stochastic approximation is a class of methods to solve intractableequations by using a sequence of approximate (and random) evaluations.See, H. Robbins and S. Monro, “A stochastic approximation method,”Annals of Mathematical Statistics, 22(3):400 407, 1951. StochasticGradient Descent (SGD) is a special type of stochastic approximationmethod that optimizes an integral using approximate gradient steps (See,Léon Bottou, “Online algorithms and stochastic approximations,” inOnline Learning and Neural Networks (David Saad, Ed., CUP 1998). It hasbeen shown that this technique is very useful in large scale learningtasks because it can provide good generalization properties with a smallnumber of passes through the data (See, Léon Bottou and OlivierBousquet, “The tradeoffs of large scale learning,” in Optimization forMachine Learning, pp. 351-368 (Suvrit Sra, Sebastian Nowozin, andStephen J. Wright, Eds, MIT Press, 2011).

The convergence properties of the SGD algorithms are directly linked tothe variance of the gradient estimate. A tradeoff between the varianceof the gradient and the convergence speed can be obtained using batching(see, for example, M. Friedlander and M. Schmidt, “Hybriddeterministic-stochastic methods for data fitting,” UBC-CS technicalreport, TR-2011-01. However, with batching, the time required for everystep increases with the size of the batch.

There remains a need for an improvement to the standard SGD algorithmsfor solving optimization problems efficiently, without introducingerrors.

BRIEF DESCRIPTION

In accordance with one aspect of the exemplary embodiment, anoptimization method includes receiving data to be sampled in optimizinga first function, which is an integral of a second function weightedaccording to a base probability distribution, the second function havingat least one first parameter to be determined in the optimizing of theintegral first function. For a number of iterations, the method includessampling the data according to a modified probability distributionfunction which depends on at least one second parameter, updating atleast one of the at least one first parameter of the second function,and updating at least one of the at least one second parameter of themodified probability distribution function. In a subsequent one of theiterations, the modified probability distribution function used insampling the data is based on the updated at least one second parameterof the distribution function. Information is output, based on the atleast one first parameter of the second function as updated in a laterone of the iterations.

In accordance with another aspect of the exemplary embodiment, a systemincludes memory which receives data to be sampled in optimizing anintegral function, which is an integral of a second function weightedaccording to a base probability distribution, the second function havingat least one first parameter to be determined in the optimizing of theintegral function. A sampling component samples the data according to amodified probability distribution function which depends on at least onesecond parameter. A first update component updates the at least onefirst parameter of the second function. A second update componentupdates the at least one second parameter of the modified probabilitydistribution function. In a subsequent one of a number of iterations,the modified probability distribution function used by the samplingcomponent for sampling the data is based on the updated second parameterof the modified probability distribution function. An information outputcomponent outputs information based on the at least one updated firstparameter of the integral function of a later one of the iterations. Aprocessor implements the first and second update components.

In accordance with another aspect of the exemplary embodiment, a methodfor optimizing a first function using adaptive weighted stochasticgradient descent is provided. The method includes receiving data to besampled in optimizing a first function which is an integral of a secondfunction depending on at least one first parameter in which elements areweighted by a probability distribution function. For a number ofiterations, the method includes sampling the data according to amodified distribution function which depends on at least one secondparameter, updating the at least one first parameter of the secondfunction based on a first step size and an estimation of the gradient ofthe second function, and updating at least one of the at least onesecond parameter of the distribution function based on a second stepsize an estimation of the variance of the gradient of the secondfunction. In a subsequent one of the iterations, the distributionfunction used in sampling the data is based on the updated at least onesecond parameter of the distribution function. Information is output,based on the at least one first parameter of the second function, asupdated in a later one of the iterations.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a functional block diagram of an exemplary computer system foroptimizing a function in accordance with one aspect of the exemplaryembodiment;

FIG. 2 is a flow chart illustrating a method for optimizing a functionin accordance with another aspect of the exemplary embodiment; and

FIG. 3 illustrates exemplary probability distribution functions.

DETAILED DESCRIPTION

The exemplary system and method provide an improvement over a standardStochastic Gradient Descent (SGD) algorithm for optimizing an integral(first) function, which is an integral of a second function of a sampledvalue and a gradient of the second function, over a space of values. Toreduce the variance of the gradient, an optimal sampling distribution islearned using a parallel SGD algorithm that learns how to minimize thevariance of the gradient of the second function.

In particular, the method accelerates an SGD algorithm for learning theparameters of the second function by learning a weight function thatminimizes the variance of the approximate gradient involved in each stepof the algorithm, while optimizing the function. While there is aconstant extra computational cost to pay at each iteration to learn theweight function, this cost is the same for each iteration, which may becompensated for by the increase of convergence speed.

With reference to FIG. 1, a functional block diagram of acomputer-implemented optimization system 10 is shown. The computersystem is configured for optimizing an integral function, given a set ofinput data 12 to be sampled, enabling parameters of the second function,and/or information based on the integral function, to be estimated. Theillustrated computer system 10 includes a processor 14, which controlsthe overall operation of the computer system 10 by execution ofprocessing instructions 16 which are stored in non-transitory memory 18,communicatively connected to the processor 14. Data memory 20 storesdata 12 being processed. Computer system 10 also includes one or moreinput/output interfaces, such as a network interface 22 and a localinterface 24. The I/O interface 24 may communicate with one or more of adisplay 26, for displaying information to users, and a user inputdevice, such as a keyboard 28 or touch or writable screen, and/or acursor control device, such as mouse, trackball, or the like, forinputting text and for communicating user input information and commandselections to the processor 12. The network interface 22 allows thecomputer to communicate input information 12 and output information 32(such as parameters of the optimized function or a solution thereof)with other devices, such as the illustrated client device 34, via acomputer network 36, such as a local area network (LAN) or wide areanetwork (WAN), or the Internet, and may comprise a modulator/demodulator(MODEM).

The various hardware components 14, 18, 20, 22, 24 of the computer 10may be all connected by a bus 38. The processor 14 executes theinstructions 16 for performing the exemplary method outlined in FIG. 2.

The computer system 10 may be a PC, such as a desktop, a laptop, palmtopcomputer, portable digital assistant (PDA), server computer, cellulartelephone, tablet computer, pager, combination thereof, or othercomputing device capable of executing instructions for performing theexemplary method.

The memory 18, 20 may represent any type of non-transitory computerreadable medium such as random access memory (RAM), read only memory(ROM), magnetic disk or tape, optical disk, flash memory, or holographicmemory. In one embodiment, the memory 18, 20 comprises a combination ofrandom access memory and read only memory. In some embodiments, theprocessor 14 and memory 18 and/or 20 may be combined in a single chip.

The digital computer processor 14 can be variously embodied, such as bya single-core processor, a dual-core processor (or more generally by amultiple-core processor), a digital processor and cooperating mathcoprocessor, a digital controller, or the like.

The exemplary system includes one or more software components 40, 42,44, 46 for implementing an algorithm for optimizing the integralfunction as described below. These components may include a samplingcomponent 40, a function parameter update component (or first updatecomponent) 42, a probability distribution function update component (orsecond update component) 44, and an information output component 46,which are best understood with reference to the method, described below.These components use the input data 12 as well as input or computed stepsizes 50 to optimize parameters 52 of the integral function 54 stored inmemory that has a set of predefined fixed parameters 55 and, in theprocess, to optimize parameters of a weighted distribution function 56used by the sampling component 40 for selecting samples of the data 12to be used in a next iteration by the function parameter updatecomponent 42.

The term “software,” as used herein, is intended to encompass anycollection or set of instructions executable by a computer or otherdigital system so as to configure the computer or other digital systemto perform the task that is the intent of the software. The term“software” as used herein is intended to encompass such instructionsstored in storage medium such as RAM, a hard disk, optical disk, or soforth, and is also intended to encompass so-called “firmware” that issoftware stored on a ROM or so forth. Such software may be organized invarious ways, and may include software components organized aslibraries, Internet-based programs stored on a remote server or soforth, source code, interpretive code, object code, directly executablecode, and so forth. It is contemplated that the software may invokesystem-level code or calls to other software residing on a server orother location to perform certain functions.

As will be appreciated, FIG. 1 is a high level functional block diagramof only a portion of the components which are incorporated into acomputer system 10. Since the configuration and operation ofprogrammable computers are well known, they will not be describedfurther.

FIG. 2 illustrates an exemplary method for optimizing an integralfunction 54, by optimizing parameters 52 for the pre-defined set offixed parameters 55, which can be performed with the system of FIG. 1.The method is best understood with reference to Algorithm 1 below. Themethod begins at S100.

At S102, a set of data points 12 to be sampled is received into memory(or simply a sample space from which the data points are to be sampled).The data points x can each be a scalar value or a vector comprising twoor more values. If not already stored, a function γ(θ) (54) to beoptimized on the data points x (12) is also received into memory. Theintegral function is defined in terms of a second function ƒ, which is afunction of x, one or more variable parameters θ (52), and a set ofconstant parameters (55), where the weight of each value x in thefunction is determined according to a base probability distribution P.

For each of a number of iterations, steps S104, S106, and S108 areperformed, e.g., in sequence. The number of iterations may be at least2, or at least 10, or at least 100, such as at least 1000.

At S104, a sampling step is performed (by the sampling component 40), bydrawing a data point x from the set X (12) of data points, where thesampling is based on a modified probability distribution (function)

which is continuous with the base probability distribution P but whichdiffers in its parameters τ (56).

At S106, parameters θ (52) of the function ƒ (54) to be optimized areupdated (by the function parameter update component 42), through astochastic gradient descent step with a magnitude specified by a stepsize parameter η_(t) (50).

At S108, parameters τ (56) of the probability distribution function

are updated (by the probability distribution function update component44), also through a stochastic gradient descent step of magnitude ε_(t)in the direction of a smaller variance for the approximate gradient. Inthe next iteration, the data point sampled in S104 is drawn based on theupdated probability distribution function

, defined by parameters τ Depending on the modification to

, this can make drawing certain samples more likely to be in thedirection of the exact gradient than in the prior iteration, whichexpedites the convergence in a manner which generally outweighs the costof the additional computation time which S108 entails.

At S110, if convergence or other stopping criterion has not beenreached, the method returns to S104 for the next iteration, otherwise toS112, where information is output (by the information output component46). The information may include one or more of the value of theoptimized function, its parameters, and information based thereon,derived from the most recent (or one of the most recent) iterations. Themethod ends at S114.

The method illustrated in FIG. 2 may be implemented in a computerprogram product that may be executed on a computer. The computer programproduct may comprise a non-transitory computer-readable recording mediumon which a control program is recorded (stored) for implementing atleast a part of the exemplary method, such as a disk, hard drive, or thelike. Common forms of non-transitory computer-readable media include,for example, floppy disks, flexible disks, hard disks, magnetic tape, orany other magnetic storage medium, CD-ROM, DVD, or any other opticalmedium, a RAM, a PROM, an EPROM, a FLASH-EPROM, or other memory chip orcartridge, or any other tangible medium from which a computer can readand use.

Alternatively, the method may be implemented in transitory media, suchas a transmittable carrier wave in which the control program is embodiedas a data signal using transmission media, such as acoustic or lightwaves, such as those generated during radio wave and infrared datacommunications, and the like.

The exemplary method may be implemented on one or more general purposecomputers, special purpose computer(s), a programmed microprocessor ormicrocontroller and peripheral integrated circuit elements, an ASIC orother integrated circuit, a digital signal processor, a hardwiredelectronic or logic circuit such as a discrete element circuit, aprogrammable logic device such as a PLD, PLA, FPGA, Graphical card CPU(GPU), or PAL, or the like. In general, any device, capable ofimplementing a finite state machine that is in turn capable ofimplementing the flowchart shown in FIG. 2, can be used to implement themethod.

As will be appreciated, the steps of the method need not all proceed inthe order illustrated and fewer, more, or different steps may beperformed.

In the following, the terms “optimization,” “minimization,” and similarphraseology are to be broadly construed as one of ordinary skill in theart would understand these terms. For example, these terms are not to beconstrued as being limited to the absolute global optimum value,absolute global minimum, and so forth. For example, minimization of afunction may employ an iterative minimization algorithm that terminatesat a stopping criterion before an absolute minimum is reached. It isalso contemplated for the optimum or minimum value to be a local optimumor local minimum value.

Stochastic Gradient Descent

First, a summary of the stochastic gradient descent method is providedby way of background.

In general, as in existing stochastic gradient methods, the exemplarymethod allows optimization, e.g., minimization, of integral functions ofthe general form:γ(θ)=∫_(X)ƒ(x;θ)dP(x)  (1)

where γ(θ) represents a possible solution,

dP(x) represents the weight applied to a sample x, where P is aprobability distribution on a measurable space Ξ from which all datasamples x are drawn, and

ƒ is a function from X×θ into

, where

is the set of real numbers and θ is a parameter from a set Θ of one ormore variable parameters (52) to be learned. The object is to identifythe values of parameters θ which provide an optimal value of γ. Avariety of different functions can be used as function ƒ, examples ofwhich are given below.

Stochastic Gradient Descent (SGD) is a stochastic approximation methodwhich involves performing approximate gradient steps equal, on averageto the true gradient ∇_(θ)γ(θ) of the function γ(θ) given distributionP, i.e., its derivative.

In some applications disclosed herein, such as logistic regressionlearning, the function ƒ is the log-likelihood and P is an empiricaldistribution

$\frac{1}{n}{\sum\limits_{i = 1}^{n}\;{\delta\left( {x,x_{i}} \right)}}$where {x₁, . . . , x_(n)} is a set of i.i.d. (independently identicallydistributed) data points sampled from an unknown distribution and δ(x,x_(i)) represents the Dirac distribution at position x_(i) (see, LéonBottou, “Online algorithms and stochastic approximations,” in OnlineLearning and Neural Networks (David Saad, Ed., CUP, 1998).

The conventional SGD method involves a sequence of iterations (steps).At a given step t, the conventional SGD method can be viewed as atwo-step procedure:

-   -   1. sample an element x_(t) in Ξ according to the distribution P;    -   2. do an approximate gradient step in which the parameter(s) are        updated: θ_(t+1)=θ_(t)−η_(t)∇_(θ)ƒ(x_(t); θ_(t)),

where η_(t) represents the gradient step size (a learning rate); and

-   -   ∇_(θ)ƒ(x_(t); θ_(t)) is the gradient of function ƒ, given data        point x_(t) and current parameters θ_(t).

One problem is in choosing the gradient step size. If η_(t) is too big,the iterations do not tend to converge. If η_(t) is too small,convergence may be very slow. In practice, when the number of uniqueelements in the support of P is small, the elements are not sampledi.i.d., but are explored one by one by doing multiple passes through thedata. For simplicity, only the random version of SGD is considered,bearing in mind that a more ordered sampling may be employed.

Weighted Stochastic Gradient Descent

Similarly to importance sampling, it can be shown that it is notnecessary to sample the elements according to the base distribution P.Rather, any distribution

>>P can be used. The notation

>>P means that P is absolutely continuous with respect to

: ∀xεX,

(x)=0→P(x)=0. This means that when

=0, P=0.

Denoting q as the density

$\frac{\mathbb{d}Q}{\mathbb{d}P},$the distribution can be written as:

$\begin{matrix}{{{\gamma(\theta)} = {\int_{X}{{\overset{\sim}{f}\left( {x;\theta} \right)}\ {\mathbb{d}{Q(x)}}}}}{where}{{\overset{\sim}{f}\left( {x;\theta} \right)} = {{{f\left( {x;\theta} \right)}\frac{\mathbb{d}{P(x)}}{\mathbb{d}{Q(x)}}} = {\frac{f\left( {x;\theta} \right)}{q(x)}.}}}} & (2)\end{matrix}$

The Weighted SGD algorithm has the following form: At each iteration:

1. sample an element x_(t) in X probabilistically, according to thedistribution

;

2. do an approximate gradient step: θ_(t+1)=θ_(t)−η_(t)∇_(θ){tilde over(ƒ)}(x_(t); θ_(t))

i.e., in the weighted SGD, the update is based on the variance of thegradient of the function ƒ.

Variance of the Approximate Gradients

By construction, unbiased estimates of the gradient ∇_(θ){tilde over(ƒ)}(x_(t); θ_(t)) can be obtained for any choice of the weightingdistribution

, i.e.,E

[∇_(θ){tilde over (ƒ)}(x _(t);θ_(t))]=∇_(θ)γ(θ).  (3)

where E

is the expectation of

.

The efficiency of the base distributions

based on their variance Σ

(θ) can be compared:

$\begin{matrix}\begin{matrix}{{\sum\limits_{Q}\;(\theta)} \equiv {{Var}_{Q}\left( {\nabla_{\theta}{\overset{\sim}{f}\left( {x_{t};\theta_{t}} \right)}} \right)}} \\{= {{E_{Q}\left\lbrack {{\nabla_{\theta}{\overset{\sim}{f}\left( {x_{t};\theta_{t}} \right)}}{\nabla_{\theta}^{T}{\overset{\sim}{f}\left( {x_{t};\theta_{t}} \right)}}} \right\rbrack} -}} \\{{E_{Q}\left\lbrack {\nabla_{\theta}{\overset{\sim}{f}\left( {x_{t};\theta_{t}} \right)}} \right\rbrack}{E_{Q}\left\lbrack {\nabla_{\theta}^{T}{\overset{\sim}{f}\left( {x_{t};\theta_{t}} \right)}} \right\rbrack}} \\{= {{E_{Q}\left\lbrack \frac{{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{\nabla_{\theta}^{T}{f\left( {x_{t};\theta_{t}} \right)}}}{{q\left( x_{t} \right)}^{2}} \right\rbrack} - {{\nabla_{\theta}{\gamma(\theta)}}{\nabla_{\theta}^{T}{\gamma(\theta)}}(5)}}}\end{matrix} & (4)\end{matrix}$

Where ∇_(θ) ^(T) . . . represents the transpose of the gradient ∇_(θ) .. . .

Now, this formula can be used to choose the optimal samplingdistribution

* that has the variance Σ

(θ) with the smallest magnitude, e.g., such that:

$\begin{matrix}{Q^{*} \in {\arg\;{\min\limits_{Q}{{tr}{\sum\limits_{Q}\;(\theta)}}}}} & (6) \\{\mspace{31mu}{= {\arg\mspace{11mu}{\min\limits_{Q}{E_{Q}\left\lbrack {\frac{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{q(x)}}^{2} \right\rbrack}}}}} & (7)\end{matrix}$

i.e., the expectation of Q is a function (norm squared) of the gradientof the variance in the function ƒ, i.e., of ∇_(θ)ƒ(x_(t); θ_(t)) anddensity q(x).

Equation (7) is obtained by noting that

$E_{Q}\left\lbrack \frac{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{q(x)} \right\rbrack$does not depend on

. Ideally it would be desirable to compute the best possible proposalfor

at each iteration, i.e., for each different value of θ_(t), i.e., toselect a distribution

on the basis of which the next sample x_(t) is drawn which minimizes thevariance in the function ƒ. However, this is not necessary, and in theexemplary method, an approximation can be performed, as follows:Adaptive Weighted Stochastic Gradient Descent

It can be demonstrated that instead of solving Equation (7) exactly, agradient step can be performed at each iteration toward lower varianceestimates.

Assume that the proposed probability distribution (from a family ofdistributions under consideration) is

_(τ) where τεT is a set of parameters for the proposed probabilitydistribution function. For example,

_(τ) could be a Gaussian distribution N(μ,Λ⁻¹) with parameters τ=(μ,Λ)ε

^(d)

S⁺, where μ is the mean and Λ is the variance of the Gaussiandistribution. S⁺ represents the set of positive definite matrices andthe symbol

means the Cartesian product.

The corresponding probability distribution function (pdf) of theproposed distribution

_(τ) at point x is denoted q(x; τ). This function is assumed to betractable, i.e. one can sample and compute q(x; τ) efficiently. Thevariance term also depends on τ:

${\Sigma_{Q}\left( {\theta,\tau} \right)} \equiv {{{Var}_{Q}\left( \frac{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{q\left( {x;\tau} \right)} \right)}.}$

An algorithm for performing stochastic gradient descent using theadaptive weighting method is shown in Algorithm 1.

Algorithm 1: Adaptive Weighted Stochastic Gradient Descent (AW-SGD)Initialize θ₀ ε Θ and τ₀ ε T. t ← 0   while not converged do (1) x_(t) □Q(. ; τ_(t))$\left. {(2)\mspace{14mu}\theta_{t + 1}}\leftarrow{\theta_{t} - {\eta_{t}\frac{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{q\left( {x_{t};\tau_{t}} \right)}}} \right.$$\left. {(3)\tau_{t + 1}}\leftarrow{\tau_{t} + {ɛ_{t}{\frac{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{q\left( {x_{t};\tau_{t}} \right)}}^{2}{\nabla_{\tau}\log}\;{q\left( {x_{t};\tau_{t}} \right)}}} \right.$(4) t ← t + 1; end while

Algorithm 1 can be summarized as follows. The algorithm is initializedwith initial values θ₀ for a set of parameters Θ of the second functionƒ to be optimized. The initial values can be selected randomly, based onexperimental data, or all set to some default value. The number t ofiterations performed is initialized to 0. Initial values for the set ofparameters τ of the modified probability distribution

are also established, e.g., by selecting mean and variance of anysuitable Gaussian distribution which is continuous with P.

In each iteration t, steps (1)-(4) of the algorithm are performed. Step(1) involves sampling a data point according to the current distribution

whose parameters τ_(t) were established in the prior iteration (S104).Eq (1) represents the sampling of a random variable according to

, which is a function of the current value of τ_(t).

Step (2) involves updating the parameters Θ of the function ƒ as afunction of the current values θ of these parameters, a gradient stepsize η_(t) and the gradient of {tilde over (ƒ)}, i.e.,

$\frac{\nabla_{\theta}{f\left( {x;\theta_{t}} \right)}}{q\left( {x;\tau_{t}} \right)}$(S106).

$\frac{\nabla_{\theta}{f\left( {x;\theta_{t}} \right)}}{q\left( {x;\tau_{t}} \right)}$is the weighted gradient of the second function, computed as a functionof θ_(t) and x_(t), where the weight is the inverse of q(x_(t); τ_(t)),so that its expectation is equal to the gradient of the second function.q(x_(t); τ_(t)) is the probability of drawing x_(t) according to themodified probability distribution function

. This step depends on τ_(t) in order to apply lower weights to highprobabilities of

in order to guarantee that the gradient in (2) is on average equal tothe gradient of the target function γ.

-   -   θ_(t) is the current value at time t of the parameter of the        second function;    -   η_(t) is the step size;    -   x_(t) is a random sample drawn according to the modified        probability distribution function; and

$\frac{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{q\left( {x_{t};\tau_{t}} \right)}.$

Step (3) involves updating the parameters τ of the distribution

as a function of a step size ε_(t) (learning rate) for the updating ofparameter τ and the gradient of the variance in Eqn. (6): i.e., thesquared norm of the ratio of the gradient of the function ƒ to q,multiplied by the gradient of the log of q:

${\frac{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{q\left( {x_{t};\tau_{t}} \right)}}^{2}{\nabla_{\tau}\log}\mspace{14mu}{{q\left( {x_{t};\tau_{t}} \right)}.}$This step aims at optimizing τ rather than θ. By changing the parametersτ, the sampling distribution Q is changed for the next iteration. As forstep (2) this is also a stochastic gradient approximation step, ratherthan an exact computation (S108).

Step (4) involves a standard step of incrementing the iteration numberto proceed to the next iteration.

Steps (1)-(4) are continued until a stopping criterion is met (S110).The iterations may be continued while γ(θ) continues to converge.However, since γ(θ) is not computed in the example algorithm, theiterations may be continued until the pre-defined fixed number ofiterations is reached. It can be shown mathematically that Algorithm 1does in fact converge, rather than achieving simply a local optimum, fora wide variety of cases. In other embodiments, the iterations arecontinued while τ continues to change by more than a threshold amount,or the like.

The step sizes η_(t) and ε_(t) can be fixed or variable (for example,decreasing as more iterations are performed). For example, η_(t) couldbe set at a fixed value of, for example, from 0.0001 to 0.2, or could beset as the minimum of a fixed value and a function of the number ofiterations, e.g.,

$\eta_{t} = {{\min\left( {\frac{1}{10},\frac{50}{t}} \right)}.}$In this example, η_(t) decreases as a function of the number ofiterations, but only begins to decrease after a predefined number (500in this example) of iterations.

Similarly, ε_(t) could be set at a fixed value of, for example, from0.0001 to 0.2, or could be set as the minimum of a fixed value and afunction of the number of iterations, e.g.,

$ɛ_{t} = {{\min\left( {\frac{1}{100},\frac{10}{t}} \right)}.}$In this example, ε_(t) decreases, but only after 1000 iterations.

It can be shown that the gradient steps in the updates are noisy (butunbiased) estimates of the gradient of the variance:

$\begin{matrix}{{\nabla_{\tau}{\Sigma_{Q}\left( {\theta,\tau} \right)}} = {\nabla_{\tau}{E_{Q}\left\lbrack {\frac{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{q\left( {x;\tau} \right)}}^{2} \right\rbrack}}} & (8) \\{\mspace{121mu}{= {E_{Q}\left\lbrack {{\frac{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{q\left( {x_{t};\tau_{t}} \right)}}^{2}{\nabla_{\tau}\log}\mspace{14mu}{q\left( {x_{t};\tau_{t}} \right)}} \right\rbrack}}} & (9)\end{matrix}$

The exemplary method avoids the need to compute Eqn. (9), by using step(3) of the algorithm.

To illustrate the sampling step (S104), FIG. 3 illustrates examples of asimplified normal base distribution P where the y axis shows theprobability of drawing a sample x from the dataset, and four differentmodified distributions

, covering the range of x values from −2 to 6. In this example, the goalis to minimize the integral of 0.5*N(x|1.0,1.0)*(x−θ)² over the realline, where N(x|m,v) is the probability density function at point x of aGaussian variable with mean in and variance v. It can be shown withstraightforward algebra that the solution can be found in closed form:the best possible θ (the one minimizing the integral) is θ=1. Inpractice, of course, the solution to many problems cannot be found inthis way. Assume that at a given step of the algorithm, θ=0.5. On eachplot, a different choice of

has been studied:

Plot 1: mean=1, variance=1;

Plot 2: mean=1, variance=1.3²;

Plot 3: mean=1.3, variance=1

Plot 4: mean=1.3, variance=1.3².

The gradient g (∇_(θ)f(x,θ)p(x)/q(x)) for each possible choice of x isalso shown for each plot. The expected gradient ∇_(θγ)(θ) is thehorizontal line at y. Formally, it is equal to the integral ofN(x|1.0,1.0)*(θ−x) over the real line. In the example, it can be shownthat it is equal on average to −0.5, as shown in the plots. However, inthe general case of stochastic gradient descent, this gradient is notknown, but uses an approximate gradient

(four different choices are shown in the figure), the estimated gradientg satisfies E [g]=0.5, where the expectation is taken with respect to

. Since the best possible choice of

is the one that give the smallest variance of the gradient, the varianceof g, Var[g] is computed and shown in the title of each plot. It can beseen that the choice q(x)=N(x|1.3,1.3) leads to a variance equal to0.77, which is smaller than the other choices of

. For this value of θ (θ=0.5), the best

distribution is not equal to the P distribution, since it gives a highervariance (Var[g]=1).

The exemplary approach is similar to adaptive importance sampling usingstochastic gradient. See, for example, W. A. Al-Qaq, M. Devetsikiotis,and J.-K. Townsend, “Stochastic gradient optimization of importancesampling for the efficient simulation of digital communication systems,”IEEE Transactions on Communications, 43(12):2975-2985, December 1995; W.A. AI-Qaq and J. K. Townsend, “A stochastic importance samplingmethodology for the efficient simulation of adaptive systems infrequency nonselective Rayleigh fading channels,” IEEE Journal onSelected Areas in Communications, 15(4):614-625, May 1997. In theseadaptive importance sampling methods, the optimal importancedistribution is learned by minimizing the variance of the samplingdistribution using a stochastic approximation procedure, but it is notused to minimize a function of the form Eqn. (1).

Adaptations to the Method

1. Batching

It is to be appreciated that Algorithm 1 can be adapted to batching(sampling multiple data points at each iteration instead of a single onein step (1)) to further reduce the variance of the gradient.

2. Scaling

To obtain good convergence properties, the gradient can be scaled,ideally in the direction of the inverse Hessian. However, this type ofsecond-order method is slow in practice. One can estimate the Hessiangreedily, as in Quasi-Newton methods such as LM-FBGS, and adapt it toSGD algorithm, similarly to the method in Friedlander.

Applications

There are many applications of SGD in solving data analytics problems.As examples, the exemplary method can be used in the following problems:

1. Large scale learning. These include classical applications of SGD.See, for example, those exemplified in Léon Bottou, “Online algorithmsand stochastic approximations,” in Online Learning and Neural Networks(David Saad, editor, Cambridge University Press, Cambridge, UK, 1998);]Léon Bottou and Olivier Bousquet, “The tradeoffs of large scalelearning,” in, Optimization for Machine Learning, pp. 351-368 (SuvritSra, Sebastian Nowozin, and Stephen J. Wright, Eds, MIT Press, 2011); M.Friedlander and M. Schmidt, “Hybrid deterministic-stochastic methods fordata fitting,” UBC-CS technical report, TR-2011-01.

2. Optimization under uncertainty (e.g., stock management).

3. Parameter estimation in latent variable models: Any application wherean expectation maximization/minimization (EM) algorithm is used. Thelikelihood can be expressed as: ∫_(Z)p(z|θ)p(x|z,θ)dz. The quantityƒ(x,θ)=p(z|θ)p(x|z,θ) is known as the complete likelihood.

4. Variational inference. SGD can be used to minimize the variationalobjective, which is sometimes intractable.

As one example, the method can be used to estimate an average futurereward, e.g., which is expressed as an integral of all future rewards.Another example is to estimate the minimum cost of stocks.

Without intending to limit the scope of the exemplary embodiment, thefollowing examples demonstrate the applicability of the method.

EXAMPLES

The exemplary algorithm can work on any integrable function, and shouldbe particularly effective on high dimensional problems, but as anillustration, the algorithm is tested on a 1D integral with oneparameter θ:γ(θ)=∫_(X) log(1+exp(θx))N(x)dx

where N (x) is the standard normal distribution (density). The integralfunction γ is convex and has therefore a unique minimum.

A conventional SGD algorithm was run using the sampling N (x) and thelearning rate optimized for this algorithm. Good performance was foundwith a learning rate

$\eta_{t} = {{\min\left( {\frac{1}{10},\frac{50}{t}} \right)}.}$Then, Algorithm 1 (adaptive weighted SGD) was run using the learningrate

$ɛ_{t} = {{\min\left( {\frac{1}{100},\frac{10}{t}} \right)}.}$

Table 1 shows numerical experiments averaged over 100 repetitions. Theerror is the distance from the optimal point computed numerically (usingquadrature). AW-SGD is the exemplary adaptive weighted SGD algorithm.SGD is the conventional stochastic gradient descent method outlinedabove (without step (3) of algorithm 1). The values in parentheses arethe standard deviations. It can be seen that the improvement is about50% of error decrease, while having an increase of only 15-30% ofcomputation time for the same number of iterations.

TABLE 1 Average time No. of (sec) to complete Method iterations errorNo. of iterations SGD 1000 0.1494 (0.096001) 0.0332 (0.013606) AW-SGD1000 0.1385 (0.095534) 0.0523 (0.009760) SGD 10,000 0.0458 (0.032266)0.2562 (0.017901) AW-SGD 10,000 0.0351 (0.024292) 0.4593 (0.017363) SGD100,000 0.0136 (0.011289) 2.5013 (0.055575) AW-SGD 100,000 0.0102(0.007616) 4.5502 (0.129989)

The standard deviation is lower with the present AW-SGD method,indicating that the method provides less variability in results, whichis advantageous.

It will be appreciated that variants of the above-disclosed and otherfeatures and functions, or alternatives thereof, may be combined intomany other different systems or applications. Various presentlyunforeseen or unanticipated alternatives, modifications, variations orimprovements therein may be subsequently made by those skilled in theart which are also intended to be encompassed by the following claims.

What is claimed is:
 1. An optimization method comprising: receiving datato be sampled in optimizing a first function, which is an integral of asecond function weighted according to a base probability distribution,the second function having at least one first parameter to be determinedin the optimizing of the integral first function; for a number ofiterations: sampling the data according to a modified probabilitydistribution function which depends on at least one second parameter;updating at least one of the at least one first parameter of the secondfunction, comprising modifying a previous value of the at least onefirst parameter of the second function as a function of a first stepsize and a weighted gradient of the second function; updating at leastone of the at least one second parameter of the modified probabilitydistribution function; wherein in a subsequent one of the iterations,the modified probability distribution function used in sampling the datais based on the updated at least one second parameter of thedistribution function; and outputting information based on the at leastone first parameter of the second function as updated in a later one ofthe iterations.
 2. The method of claim 1, wherein the updating of the atleast one first parameter of the second function is performed bygradient descent.
 3. The method of claim 1, wherein the updating of theat least one second parameter of the modified probability distributionfunction is performed by gradient descent.
 4. The method of claim 1,wherein the modified probability distribution function comprises anormal distribution function and the at least one second parameter whichis updated comprises at least one of its mean and variance.
 5. Themethod of claim 1, wherein the first step size has a value which is afunction of a number of the iterations which have been performed.
 6. Themethod of claim 1, wherein the updating of the at least one firstparameter of the second function is performed according to the equation:$\left. \theta_{t + 1}\leftarrow{\theta_{t} - {\eta_{t}\frac{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{q\left( {x_{t};\tau_{t}} \right)}}} \right.,$where: θ_(t) is the current value at iteration t of the parameter of thesecond function; η_(t) is the step size; x_(t) is the value of a randomsample drawn from the data according to the modified probabilitydistribution function; and$\frac{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{q\left( {x_{t};\tau_{t}} \right)}$is the weighted gradient of the second function, computed as a functionof θ_(t) and x_(t), weighted by the inverse of q(x_(t);τ_(t)), and whereq(x_(t);τ_(t)) is the probability of drawing x_(t) according to themodified probability distribution function.
 7. The method of claim 1,wherein the updating of the at least one second parameter of themodified probability distribution function minimizes the variance of theapproximate gradient in the second function in each iteration.
 8. Themethod of claim 1, wherein the integral of the function to be optimizedis convex.
 9. The method of claim 1, wherein the modified probabilitydistribution function differs from the base probability distribution andis absolutely continuous with respect to the base distribution function.10. The method of claim 1, wherein the updating of the at least onefirst parameter of the second function is based on a first step size andan estimation of the gradient of the second function and the updating ofthe at least one of the at least one second parameter of thedistribution function based on a second step size an estimation of thevariance of the gradient of the second function.
 11. The method of claim1, wherein the output information comprises at least one of: a solutionto the integral function; the updated at least one parameter of thesecond function; and information based thereon.
 12. An optimizationmethod comprising: receiving data to be sampled in optimizing a firstfunction, which is an integral of a second function weighted accordingto a base probability distribution, the second function having at leastone first parameter to be determined in the optimizing of the integralfirst function; for a number of iterations: sampling the data accordingto a modified probability distribution function which depends on atleast one second parameter; updating at least one of the at least onefirst parameter of the second function; updating at least one of the atleast one second parameter of the modified probability distributionfunction, comprising modifying a previous value of the at least onesecond parameter as a function of a second step size and a quantityequal on average to a gradient of the variance of the first parameterupdate; wherein in a subsequent one of the iterations, the modifiedprobability distribution function used in sampling the data is based onthe updated at least one second parameter of the distribution function;and outputting information based on the at least one first parameter ofthe second function as updated in a later one of the iterations.
 13. Themethod of claim 12, wherein the updating of the at least one firstparameter of the second function comprises modifying a previous value ofthe at least one first parameter of the second function as a function ofa first step size and a weighted gradient of the second function. 14.The method of claim 12, wherein the second step size has a value whichis a function of a number of the iterations which have been performed.15. The method of claim 12, wherein the updating of the at least one ofthe at least one second parameter is performed according to theequation:$\left. \tau_{t + 1}\leftarrow{\tau_{t} + {ɛ_{t}{\frac{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{q\left( {x_{t};\tau_{t}} \right)}}^{2}{\nabla_{\tau}\log}\mspace{14mu}{q\left( {x_{t};\tau_{t}} \right)}}} \right.,$where: τ_(t) is the current value of the second parameter of themodified probability distribution function at iteration t, θ_(t) is thecurrent value of the first parameter of the second function; ε_(t) isthe step size; x_(t) is the value of a random sample drawn from the dataaccording to the modified probability distribution function; and${\frac{\nabla_{\theta}{f\left( {x_{t};\theta_{t}} \right)}}{q\left( {x_{t};\tau_{t}} \right)}}^{2}{\nabla_{\tau}\log}\mspace{14mu}{q\left( {x_{t};\tau_{t}} \right)}$is equal on average to the gradient of the variance of the firstparameter update.
 16. An optimization method comprising: receiving datato be sampled in optimizing a first function, which is an integral of asecond function weighted according to a base probability distribution,the second function having at least one first parameter to be determinedin the optimizing of the integral first function, the integral functionbeing of the form:γ(θ)=∫_(X)θ(x;θ)dP(x)  (1), where: dP(x) represents the weight appliedto a sample x, and where P is a probability distribution on a measurablespace Ξ from which the data samples are drawn, ƒ represents the secondfunction; and θ represents the at least one first parameter; for anumber of iterations: sampling the data according to a modifiedprobability distribution function which depends on at least one secondparameter; updating at least one of the at least one first parameter ofthe second function; updating at least one of the at least one secondparameter of the modified probability distribution function; wherein ina subsequent one of the iterations, the modified probabilitydistribution function used in sampling the data is based on the updatedat least one second parameter of the distribution function; andoutputting information based on the at least one first parameter of thesecond function as updated in a later one of the iterations.
 17. Acomputer program product comprising a non-transitory recording mediumstoring instructions, which when executed on a computer causes thecomputer to perform the method of claim
 1. 18. A system comprisingmemory which stores instructions for performing the method of claim 1and a processor, in communication with the memory, which executes theinstructions.
 19. A system comprising: memory which receives data to besampled in optimizing an integral function, which is an integral of asecond function weighted according to a base probability distribution,the second function having at least one first parameter to be determinedin the optimizing of the integral function; a sampling component forsampling the data according to a modified probability distributionfunction which depends on at least one second parameter; a first updatecomponent for updating the at least one first parameter of the secondfunction; a second update component for updating the at least one secondparameter of the modified probability distribution function; wherein atleast one of: the updating of the at least one first parameter of thesecond function is based on a first step size and an estimation of thegradient of the second function, and the updating of the at least one ofthe at least one second parameter of the distribution function based ona second step size an estimation of the variance of the gradient of thesecond function; wherein in a subsequent one of a number of iterations,the modified probability distribution function used by the samplingcomponent for sampling the data is based on the updated second parameterof the modified probability distribution function; an information outputcomponent for outputting information based on the at least one updatedfirst parameter of the integral function of a later one of theiterations; and a processor which implements the first and second updatecomponents.
 20. A method for optimizing a function using adaptiveweighted stochastic gradient descent, comprising: receiving data to besampled in optimizing a first function which is an integral of a secondfunction depending on at least one first parameter in which elements areweighted by a base probability distribution function; for a number ofiterations: sampling the data according to a modified probabilitydistribution function which depends on at least one second parameter,wherein the modified probability distribution function differs from thebase probability distribution function and is absolutely continuous withrespect to the base probability distribution function; updating the atleast one first parameter of the second function based on a first stepsize and an estimation of the gradient of the second function; updatingat least one of the at least one second parameter of the distributionfunction based on a second step size an estimation of the variance ofthe gradient of the second function; wherein in a subsequent one of theiterations, the distribution function used in sampling the data is basedon the updated at least one second parameter of the distributionfunction; outputting information based on the at least one firstparameter of the second function as updated in a later one of theiterations.